The inner product, also known as the dot product can be computed in R using the following function:

Given the two vectors:

u <- c(2, -5, -1)
v <- c(3, 2, -3)

t(u) %*% v
##      [,1]
## [1,]   -1
dot_product <- function(vector1, vector2) {
  t(vector1) %*% vector2
}

dot_product(u, v)
##      [,1]
## [1,]   -1

The length (or norm) of a vector is the square root the sum of all entries in the vector squared.

A vector whose length is 1 is called a unit vector. Example:

Let v = (1, -2, 2, 0). Find a unit vector u in the same direction as v.

Compute the length of v:

v <- c(1, -2, 2, 0)

length_vector <- function(vector, norm = FALSE) {

length_v <- sqrt(sum(vector*vector))

## Multiply __v__ by 1/length_v

  if (norm){
  
  return((1/length_v)*vector)
  
  }
  
  else {
    return(length_v)
  }
}

#Length of a vector:
length_vector(v)
## [1] 3
#Unit Vector
length_vector(v, TRUE)
## [1]  0.3333333 -0.6666667  0.6666667  0.0000000
# Check that vector is a unit vector
# It returns 1 so this is true. 
length_vector(length_vector(v, TRUE))
## [1] 1

For u and v in Rn, the distance between u and v, written as dist(u, v), is the length of the vector u-v. That is, dist(u, v) = ||u - v||

Example: Compue the distance between the vectors u = (7,1) and v = (3,2)

# Let 
u <- c(7, 1)
v <- c(3, 2)
u_minus_v <- u - v
length_vector(u_minus_v)
## [1] 4.123106

Two vectors are orthogonal to each other if the dot product is 0.

u <- c(0, 1)
v <- c(1, 0)
dot_product(u, v)
##      [,1]
## [1,]    0

Orthogonal Sets

u1 <- c(3, 1, 1)
u2 <- c(-1, 2, 1)
u3 <- c(-1/2, -2, 7/2)

dot_product(u1, u2)
##      [,1]
## [1,]    0
dot_product(u2, u3)
##      [,1]
## [1,]    0
dot_product(u1, u3)
##      [,1]
## [1,]    0

An orthogonal basis for a subsace W of Rn is a basis for W that is also an orthogonal set.

y <- c(6, 1, -8)

# Returns y
dot_product(y, u1)/dot_product(u1, u1)*u1+dot_product(y, u2)/dot_product(u2, u2)*u2 + dot_product(y, u3)/dot_product(u3, u3)*u3
## [1]  6  1 -8

Orthogonal projections

y <- c(7, 6)
u <- c(4, 2)

orthogonal_proj <- function(y, Un, y_minus_yhat = FALSE) {
  y_hat <- dot_product(y,Un)/dot_product(Un,Un)*Un
  if(y_minus_yhat) {
    return(y - y_hat)
  }
  else {
    return(y_hat)
  }
  
}

# Orthogonal projection of y onto u
orthogonal_proj(y, u)
## [1] 8 4
# Component of y orthogonal to u
orthogonal_proj(y, u, TRUE)
## [1] -1  2
# Verify that they sum to y
orthogonal_proj(y, u) + orthogonal_proj(y, u, TRUE)
## [1] 7 6

# Verify that they are perpindicular
dot_product(orthogonal_proj(y, u), orthogonal_proj(y, u, TRUE))
##      [,1]
## [1,]    0

length_vector(orthogonal_proj(y, u, TRUE))
## [1] 2.236068

v1 <- c(3/sqrt(11), 1/sqrt(11), 1/sqrt(11))
v2 <- c(-1/sqrt(6), 2/sqrt(6), 1/sqrt(6))
v3 <- c(-1/sqrt(66), -4/sqrt(66), 7/sqrt(66))

dot_product(v1, v2)
##      [,1]
## [1,]    0
dot_product(v1, v3)
##              [,1]
## [1,] 5.551115e-17
dot_product(v2, v3)
##      [,1]
## [1,]    0
dot_product(v1, v1)
##      [,1]
## [1,]    1
dot_product(v2, v2)
##      [,1]
## [1,]    1
dot_product(v3, v3)
##      [,1]
## [1,]    1

U <- matrix(c(v1, v2, v3), ncol = 3)
round(t(U)%*%U, 1)
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1
# I used round to show that this is nearly the identity matrix.
t(U)%*%U
##              [,1]         [,2]         [,3]
## [1,] 1.000000e+00 7.681559e-19 3.007996e-17
## [2,] 7.681559e-19 1.000000e+00 2.696773e-17
## [3,] 3.007996e-17 2.696773e-17 1.000000e+00

u1 <- c(2, 5, -1)
u2 <- c(-2, 1, 1)
y <- c(1, 2, 3)

# The orthogonal projection of y onto W
orthogonal_proj(y, u1) + orthogonal_proj(y, u2)
## [1] -0.4  2.0  0.2
y - (orthogonal_proj(y, u1) + orthogonal_proj(y, u2))
## [1] 1.4 0.0 2.8

u1 <- c(2, 5, -1)
u2 <- c(-2, 1, 1)
y <- c(1, 2, 3)

# W = Span {u1, u2}

# The closest point in W to y is

orthogonal_proj(y, u1) + orthogonal_proj(y, u2)
## [1] -0.4  2.0  0.2

y <- c(-1, -5, 10)
u1 <- c(5, -2, 1)
u2 <- c(1, 2, -1)

# The distance from y to W is ||y-y_hat||

y - orthogonal_proj(y, u1) - orthogonal_proj(y, u2)
## [1] 0 3 6
length_vector(y - orthogonal_proj(y, u1) - orthogonal_proj(y, u2))
## [1] 6.708204

[](

library(pracma)

x1 <- c(1, 1, 1, 1)
x2 <- c(0, 1, 1, 1)
x3 <- c(0, 0, 1, 1)

# Construct an orthogonal basis for W

# Step 1: 
v1 <- x1

# Step 2:
v2 <- x2 - orthogonal_proj(x2, v1)
v1
## [1] 1 1 1 1
v2
## [1] -0.75  0.25  0.25  0.25
#Step 3

v3 <- x3 - orthogonal_proj(x3, v1) - orthogonal_proj(x3, v2)
v3
## [1]  0.0000000 -0.6666667  0.3333333  0.3333333
A <- matrix(c(x1, x2, x3), ncol = 3)
gS <- gramSchmidt(A)

round(gS$Q %*% gS$R, 1)
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    1    1    0
## [3,]    1    1    1
## [4,]    1    1    1

v1 <- c(3, 6, 0)
v2 <- c(0, 0, 2)

u1 <- v1/length_vector(v1)
u2 <- v2/length_vector(v2)
u1
## [1] 0.4472136 0.8944272 0.0000000
u2
## [1] 0 0 1

v1 <- c(1, 1, 1, 1)
v2 <- c(-3, 1, 1, 1)
v3 <- c(0, -2/3, 1/3, 1/3)

u1 <- v1/length_vector(v1)
u2 <- v2/length_vector(v2)
u3 <- v3/length_vector(v2)

Q <- matrix(c(u1, u2, u3), ncol = 3)

Q
##      [,1]       [,2]        [,3]
## [1,]  0.5 -0.8660254  0.00000000
## [2,]  0.5  0.2886751 -0.19245009
## [3,]  0.5  0.2886751  0.09622504
## [4,]  0.5  0.2886751  0.09622504
A
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    1    1    0
## [3,]    1    1    1
## [4,]    1    1    1
R <- t(Q) %*% A
R
##              [,1]      [,2]      [,3]
## [1,] 2.000000e+00 1.5000000 1.0000000
## [2,] 1.110223e-16 0.8660254 0.5773503
## [3,] 0.000000e+00 0.0000000 0.1924501
round(R, 7)
##      [,1]      [,2]      [,3]
## [1,]    2 1.5000000 1.0000000
## [2,]    0 0.8660254 0.5773503
## [3,]    0 0.0000000 0.1924501
gramSchmidt(A)
## $Q
##      [,1]       [,2]       [,3]
## [1,]  0.5 -0.8660254  0.0000000
## [2,]  0.5  0.2886751 -0.8164966
## [3,]  0.5  0.2886751  0.4082483
## [4,]  0.5  0.2886751  0.4082483
## 
## $R
##      [,1]      [,2]      [,3]
## [1,]    2 1.5000000 1.0000000
## [2,]    0 0.8660254 0.5773503
## [3,]    0 0.0000000 0.8164966

A <- matrix(c(4, 0, 0, 2, 1, 1), ncol = 2, byrow = TRUE)

b <- c(2, 0, 11)

solve(t(A)%*%A)%*%t(A)%*%b
##      [,1]
## [1,]    1
## [2,]    2
least_squares_sol <- function(A, b) {
  solve(t(A)%*%A)%*%t(A)%*%b
}
least_squares_sol(A, b)
##      [,1]
## [1,]    1
## [2,]    2

# Least squares error:

b <- c(2, 0, 11)
A <- matrix(c(4, 0, 0, 2, 1, 1), ncol = 2, byrow = TRUE)
x_hat <- c(1, 2)

length_vector(b - A%*%x_hat)
## [1] 9.165151

A <- matrix(c(1, 3, 5, 1, 1, 0, 1, 1, 2, 1, 3, 3), ncol = 3, byrow = TRUE)
b <- c(3, 5, 7, -3)

least_squares_sol_gs <- function(A, b) {

gS <- gramSchmidt(A)

solve(gS$R)%*%t(gS$Q)%*%b
}

least_squares_sol_gs(A, b)
##      [,1]
## [1,]   10
## [2,]   -6
## [3,]    2
least_squares_sol(A, b)
##      [,1]
## [1,]   10
## [2,]   -6
## [3,]    2

# Fit a line

# Two arbitrary vectors of equal length
x <- c(0, 1, 2, 3, 4, 5, -2, 3, 5)
y <- c(0, 1, 4, 9, 16, 25, 0, 3, 7)


# Matrix Computation
A <- matrix(c(rep(1, length(x)),x), ncol = 2)
A
##       [,1] [,2]
##  [1,]    1    0
##  [2,]    1    1
##  [3,]    1    2
##  [4,]    1    3
##  [5,]    1    4
##  [6,]    1    5
##  [7,]    1   -2
##  [8,]    1    3
##  [9,]    1    5
AtA <- t(A)%*%A
B <- t(A) %*% y
vector <- solve(AtA)%*%B
vector
##          [,1]
## [1,] 1.000000
## [2,] 2.666667
# Plot points and line
slope <- vector[2,]
intercept <- vector[1,]

plot(x, y)
curve(intercept+slope*x, add = TRUE)

# Fit a parabola

# Two arbitrary vectors of equal length
x <- c(0, 1, 2, 3, 4, 5, -2, 3, 5)
y <- c(0, 1, 4, 9, 16, 25, 0, 3, 7)


# Matrix Computation
A <- matrix(c(rep(1, length(x)),x, x^2), ncol = 3)
A
##       [,1] [,2] [,3]
##  [1,]    1    0    0
##  [2,]    1    1    1
##  [3,]    1    2    4
##  [4,]    1    3    9
##  [5,]    1    4   16
##  [6,]    1    5   25
##  [7,]    1   -2    4
##  [8,]    1    3    9
##  [9,]    1    5   25
AtA <- t(A)%*%A
B <- t(A) %*% y
vector <- solve(AtA)%*%B
vector
##            [,1]
## [1,] -0.0998308
## [2,]  0.9949239
## [3,]  0.4839255
# Plot points and line
slope <- vector[2,]
intercept <- vector[1,]
squared <- vector[3,]
plot(x, y)
curve(intercept+slope*x+squared*x*x, add = TRUE)

# Fit anything?

# Two arbitrary vectors of equal length
x <- c(0, 1, 2, 3, 4, 5, -2, 3, 5)
y <- c(0, 1, 4, 9, 16, 25, 0, 3, 7)


# Matrix Computation
A <- matrix(c(rep(1, length(x)),x, x^2), ncol = 3)
A
##       [,1] [,2] [,3]
##  [1,]    1    0    0
##  [2,]    1    1    1
##  [3,]    1    2    4
##  [4,]    1    3    9
##  [5,]    1    4   16
##  [6,]    1    5   25
##  [7,]    1   -2    4
##  [8,]    1    3    9
##  [9,]    1    5   25
AtA <- t(A)%*%A
B <- t(A) %*% y
vector <- solve(AtA)%*%B
vector
##            [,1]
## [1,] -0.0998308
## [2,]  0.9949239
## [3,]  0.4839255
# Plot points and line
slope <- vector[2,]
intercept <- vector[1,]
squared <- vector[3,]
plot(x, y)
curve(intercept+slope*x+squared*x*x, add = TRUE)

X <- matrix(c(1, 2, 1, 5, 1, 7, 1, 8), ncol = 2, byrow = TRUE)
y <- c(1, 2, 3, 3)

least_squares_line <- function(X, y) {
coefficients <- solve(t(X)%*%X)%*%t(X)%*%y

plot(X[,2], y)  
curve(coefficients[1,] + coefficients[2,]*x, add = TRUE)
}

least_squares_line(X, y)

#Playing with some dui data

dui_month <- read.csv("C:/Users/Freddy/Documents/dui_month.csv", stringsAsFactors = FALSE, header = TRUE)
dui_month <- dui_month[,-1]
dui_month <- t(dui_month)
dui_month
##     [,1] [,2] [,3] [,4] [,5]
## Jan  464  459  548  502  555
## Feb  421  431  559  514  509
## Mar  455  573  609  592  484
## Apr  412  388  519  522  481
## May  418  495  539  553  445
## Jun  352  527  533  484  474
## Jul  395  489  512  486  479
## Aug  393  547  547  539  542
## Sep  449  500  491  441  484
## Oct  451  511  499  539  525
## Nov  368  527  488  492  481
## Dec  487  560  525  491  385
length(dui_month)
## [1] 60
plot(1:60, dui_month)

dui_month
##     [,1] [,2] [,3] [,4] [,5]
## Jan  464  459  548  502  555
## Feb  421  431  559  514  509
## Mar  455  573  609  592  484
## Apr  412  388  519  522  481
## May  418  495  539  553  445
## Jun  352  527  533  484  474
## Jul  395  489  512  486  479
## Aug  393  547  547  539  542
## Sep  449  500  491  441  484
## Oct  451  511  499  539  525
## Nov  368  527  488  492  481
## Dec  487  560  525  491  385
dui <- matrix(dui_month, ncol = 1, byrow = FALSE)
plot(1:60, dui)

X <- matrix(c(rep(1, 60), 1:60, (1:60)^2, (1:60)^3), ncol = 4)

least_squares_line <- function(X, y) {
  coefficients <- solve(t(X)%*%X)%*%t(X)%*%y
  
  plot(X[,2], y)  
  curve(coefficients[1,] + coefficients[2,]*x + coefficients[3,]*x^2 + coefficients[4,]*x^3, add = TRUE)
}

least_squares_line(X, dui)